/*
 * main.cpp
 *
 *  Created on: Jun 18, 2011
 *      Author: blyth
 */

#include <DataWriter.h>
#include <SimDataReader.h>
#include <PeakReference.h>
#include <crystalMap.h>
#include <iostream>
#include <globals.hh>
#include "CLHEP/Random/MTwistEngine.h"
#include "CLHEP/Random/RandGauss.h"
using namespace std;

int mround(double in) {
	int whole = in;
	double remain = in - whole;
	if (remain >= 0.5)
		return whole + 1;
	return whole;
}

int main(int argc, char** argv) {
	if (argc != 3) {
		cout << "Please give simulated data input as first argument, and run number as the second" << endl;
		return 0;
	}

	unsigned int const sourceRefSize = 44;
	double* sourceRef =
			new double[sourceRefSize] { 0.00229951, 0.0217842, 0.0613074, 0.156733, 0.305301, 0.49182, 0.68358,
					0.845512, 0.950986, 0.984521, 0.944636, 0.842791, 0.687357, 0.527077, 0.380956, 0.260087, 0.168696,
					0.105998, 0.0685462, 0.0484331, 0.0353872, 0.0332753, 0.0334183, 0.0331892, 0.0318952, 0.0299273,
					0.0283969, 0.0277491, 0.0273934, 0.0265412, 0.0250565, 0.0232322, 0.0213166, 0.0190748, 0.0159684,
					0.0131479, 0.0104398, 0.00725754, 0.00329131, 0.00210918, 0.00147583, 0.000903992, 0.000475297,
					0.000122942 };
	for (int i = 0; i < sourceRefSize; i++)
		sourceRef[i] *= 0.111241333;

	PeakReference pr(sourceRef, sourceRefSize);
	TH2D* dataHisto = new TH2D("dataHisto", "dataHisto", 72, -0.9, 0.9, 72, -0.9, 0.9);
	SimDataReader sdr(argv[1]);
	DataWriter dw(atoi(argv[2]));

	CLHEP::MTwistEngine randEngine(43867547654);
	CLHEP::RandGauss randGauss(randEngine);

	unsigned int** dout = new unsigned int*[2760];
	double* sout = new double[2760];
	for (int c = 0; c < 2760; c++)
		dout[c] = new unsigned int[49];

	//	int count = 0;
	while (sdr.ReadNext(dataHisto)/* && count < 1000*/) {
		for (int c = 0; c < 2760; c++) {
			double x = crystalMap[c][0];
			double y = crystalMap[c][1];
			double sum;
			int i = x / 25 + 36.5;
			int j = y / 25 + 36.5;
			if (c < 2240) {
				sum = dataHisto->GetBinContent(i, j) / MeV;
			} else {
				sum = dataHisto->GetBinContent(i, j) / MeV;
				sum += dataHisto->GetBinContent(i + 1, j + 1) / MeV;
				sum += dataHisto->GetBinContent(i + 1, j) / MeV;
				sum += dataHisto->GetBinContent(i, j + 1) / MeV;
			}
			sout[c] = sum;
			dout[c][48] = randEngine.flat() * 1008 + 16;
			for (int k = 0; k < 16; k++) {
				int roundedNoise = mround((int) randGauss.fire(0., 2.65));
				if (roundedNoise < -16)
					roundedNoise = -16;
				else if (roundedNoise > 16)
					roundedNoise = 16;
				dout[c][k] = dout[c][48] + roundedNoise;
			}
			double shift = randEngine.flat() * 0.5;
			for (int k = 16; k < 48; k++) {
				int roundedNoise = mround((int) randGauss.fire(0., 2.65));
				if (roundedNoise < -16)
					roundedNoise = -16;
				else if (roundedNoise > 16)
					roundedNoise = 16;
				dout[c][k] = dout[c][48] + roundedNoise + mround(
						sum * 155.2236 * pr.cubicInterpolationAt(k - 16 + shift));
				if (dout[c][k] > 16384)
					dout[c][k] = 16384;
			}
		}
		dw.writeData(dout, sout);
		delete dataHisto;
		dataHisto = new TH2D("dataHisto", "dataHisto", 72, -0.9, 0.9, 72, -0.9, 0.9);
		//	count++;
	}

	delete[] sout;
	for (int c = 0; c < 2760; c++)
		delete[] dout[c];
	delete[] dout;
	delete dataHisto;
	delete[] sourceRef;

	return 0;
}
